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1 Summary 


A numerical ocean tide model has been developed and tested using highly accurate 
TOPEX/Poseidon (T/P) tidal solutions. The hydrodynamic model is based on time 
stepping a finite difference approximation to the non-linear shallow water equations. Two 
novel features of our implementation are a rigorous treatment of self attraction and load- 
ing (SAL), and a physically based parameterization for internal tide (IT) radiation drag. 
The model was run for a range of grid resolutions, and with variations in model parame- 
ters and bathymetry. For a rational treatment of SAL and IT drag, the model run at high 
resolution (1/12 degree) fits the T/P solutions to within 5 cm RMS in the open ocean. 
Both the rigorous SAL treatment and the IT drag parameterization are required to obtain 
solutions of this quality. The sensitivity of the solution to perturbations in bathymetry 
suggest that the fit to T/P is probably now limited by errors in this critical input. Since 
the model is not constrained by any data, we can test the effect of dropping sea-level to 
match estimated bathymetry from the last glacial maximum (LGM). Our results suggest 
that the « 100 m drop in sea-level in the LGM would have significantly increased tidal 
amplitudes in the North Atlantic, and increased overall tidal dissipation by about 40%. 
However, details in tidal solutions for the past 20 ka are sensitive to the assumed strati- 
fication. IT drag accounts for a significant fraction of dissipation, especially in the LGM 
when large areas of present day shallow sea were exposed, and this parameter is poorly 
constrained at present. 


2 Introduction 


The evolution of the mooon’s orbit around the Earth is closely linked to the history of 
tidal dissipation in the ocean (e.g., Bills, and Ray, 1999). This connection, and the clear 
evidence that tidal dissipation must have varied significantly over time, have spurred a 
number of efforts to model the ocean tides in the recent and distant past (e.g., Thomas 
and Sundermann; 1999; and references therein). It is far from clear how reliable these 
modeling studies are, given that bathymetry in the past is only known approximately, 
and given the challenges in numerical modeling of even the modern day tides. 

With the availability of TOPEX/Poseidon (T/P) altimeter data, barotropic tidal eleva- 
tions in the open ocean are now known to within roughly 1 cm (Shum et al., 2000). 
However, to obtain solutions with this accuracy it has been necessary to use some sort of 
data assimilation or empirical mapping approach (e.g., Egbert et ah, 1994; Schrama and 
B'&yj 1994). The results of Le Provost et ah (1994) with the finite element hydrodynamic 
model FES94 suggested that accuracies approaching the empirically constrained solutions 
could be obtained with a purely hydrodynamic model, if high enough numerical resolution 
were used. However, in this work modeling was accomplished separately for individual 
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ocean basins, and the results were spliced together. Open boundary conditions for each 
basin were adjusted to achieve agreement with tide gauges, sometimes even using a strong 
constraint data assimilation approach (Lyard and Genco, 1994), so these solutions were 
not really independent of observational data. Subsequent improvements to the Grenoble 
modeling code have allowed simultaneous modeling of the full globe with no adjustable 
open boundary conditions. The resulting global solutions fit validation tide gauges much 
more poorly (F. Lyard, personal communication), demonstrating the critical role played 
by data in FES94. 

There is now strong evidence that significant energy is transfered from barotropic to 
internal tides over rough topography in the deep ocean (Egbert and Ray, 2000; 2001). 
Since this “internal tide drag” appears to account for about 1/3 of the total tidal dissi- 
pation, accurate modeling of the barotropic tide almost certainly will require accounting 
for this conversion, either by modeling the full three-dimensional stratified ocean, or by 
some sort of parameterization (Jayne and St. Laurent, 2001). Only the later approach is 
computationally feasible at this time. 

In this report we describe our efforts to develop a global barotropic hydrodynamic model 
which reproduces the present day tidal elevation fields using only the astronomical forcing, 
with no data assimilation or boundary condition constraints. We use a straightforward 
modeling approach, based on finite difference time-stepping of the non-linear barotropic 
shallow water equations (SWE), with a rigorous treatment of ocean self attraction and 
loading (SAL) effects, and several parameterizations of internal tide (IT) drag. Computa- 
tions are done for a wide range of nearly global model grids, with resolutions ranging from 
1 -1/12 . We also do experiments with small perturbations to our standard bathymetry, 
to test sensitivity of solutions to this critical input parameter. For all of the results re- 
ported here we focus on the principal lunar constituent M 2 , but other constituents were 
included in the modeling effort. As we shall show, including a parameterization for IT 
drag significantly improves the fidelity of the solution. The best model results (5 cm 
RMS misfit between deep-ocean T/P elevations and the numerical model) are obtained 
with the highest resolution grid, and with rigorous treatment of ocean self-attraction and 
loading. The level of accuracy achieved is in fact reasonably consistent with the effect of 
likely errors in the currently best available bathymetry. 

Finally, we consider the effect of dropping sea-level to that estimated for the last glacial 
maximum (LGM). The IT drag parameterization we use depends on ocean stratification, 
and it is not clear how this effect should be modeled for the LGM. We test several 
scenarios, including stratification similar to the modern ocean, and significant reductions 
and increases in stratification. As we shall show, the effects of LGM sea-level changes 
are much larger than errors in the modern solutions. In particular, tidal dissipation 
increases by up to 50% in the LGM. However, the sensitivity of the model results to both 
bathymetry and stratification suggests that efforts to accurately model tidal dissipation 
in even the recent past remain very challenging. 
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3 Hydrodynamic Modeling 

3.1 Finite-difference SWE 


We assume shallow water dynamics of the form 


<9U 

dt 


+ f xU + U- VU + a ff V 2 U + ^V(C-Cs AL )+J c ' = fo 


( 1 ) 


<K 

dt 


-v-u. 


(2) 


Here is the tidal elevation; U is the volume transport vector, equal to velocity times 
water depth H; f is the Coriolis parameter (oriented to the local vertical), T is the fric- 
tional or dissipative stress, and the term a#V 2 U is a crude parameterization of horizontal 
turbulent eddy viscosity, included primarily to improve numerical stability. For most of 
the results discussed here we have taken a H = 10 3 ra 2 s _1 . Elevations in the open ocean 
were found to be insensitive to the exact value of this parameter unless it was signifi- 
cantly increased (e.g., to 10 3 m 2 s ' ) . The global solution is also only weakly dependent 
on inclusion of the non-linear terms in (1). The astronomical tide generating force, which 
includes the Earth’s body tide (Hendershott, 1972) is denoted by f 0 . We included up 
to 8 constituents (M 2 , S 2 , N 2 , K 2 , Ki, Oi, Pi Q\) in fo. However, very similar results 
were obtained (for M 2 ) when the forcing was restricted to the dominant semi-diurnal and 
diurnal constituents M 2 and Ki, so only these two constituents were used for most of the 
extensive modeling experiments described here. 


Tidal loading and self-attraction (Hendershott, 1972; Ray, 1998) are accounted for by the 
term Csal, which we discuss in detail below. We solve the system of equations (1) and (2) 
numerically on a C-grid, following the finite difference time-stepping approach described 
in Egbert et al. (1994). 

All grids were nearly global, covering the area from 86°S to 82.25°N. Bathymetry was 
interpolated (and smoothed where appropriate) from a combination of the 1/30 degree 
database of Smith and Sandwell (1997) in deep water equatorward of 72°, ET0P05 
(National Geophysical Data Center, 1992) in shallow areas and the Arctic, and a new 
compilation of bathymetry for the Antarctic (L. Padman, personal communication). Open 
boundary conditions at the top of the domain in the Arctic were specified elevations, taken 
from the global FES94 solution (Le Provost et al., 1994), Tests with variants (including 
the Arctic assimilation solutions of Kivman (1997)) showed that elevations outside the 
arctic were nearly independent of the details in these boundary conditions. Boundary 
conditions at the coast were zero normal flow, and no slip. 

The dissipation term T — Tb + Pit included components for bottom boundary layer 
drag and IT wave radiation. The first (and standard) component was parameterized in 
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the usual way as quadratic in velocity 


T B = (c D \\v\\/H)\J, (3) 

where v is the total velocity vector (in particular, including all tidal constituents), and 
the value of the non-dimensional parameter c D is approximately 0.0025. Test with a range 
of values for c D showed that deep-water elevations were only weakly sensitive to the value 
of this parameter. Several different parameterizations of IT wave drag were tried. 


3.2 Parameterization of Internal Tide Drag 


Jayne and St. Laurent (2001) describe a simple parameterization of IT drag based on 
a linear analysis of energy flux into the internal wave field due originally to Bell (1975). 
In this theory energy conversion from barotropic to baroclinic tides by small amplitude 
sinusoidal topography of amplitude h and wavenumber k is 


( w 2 _ f 2U/2 

Ef = — p 0 Kh 2 Nu 2 , (4) 

where N is the buoyancy frequency, p 0 the mean ocean density, / the Coriolis parameter, u 
the tidal frequency, and u the barotropic tidal velocity perpendicular to the topography. 
Since the theory is linear, superposition may be used to estimate conversion for more 
complex topographic variations. To apply (4) Jayne and St. Laurent (2001) estimated 
the height of the scattering topography as the RMS of bathymetric variations not resolved 
by their 1/2 numerical grid, and, ignoring the dependence on tidal frequency ui, obtained 
a spatially varying linear drag coefficient 


Cit 


1 

2 H 


nh 2 N. 


( 5 ) 


The buoyancy frequency N (for the ocean bottom) was obtained from Levitus (1999), the 
wavenumber k was left as a tunable parameter, and the linear dissipation term Tit = c IT U 
was then added to the quadratic bottom drag (3). Note that c IT is just the linear drag 
coefficient required to match the energy loss to the barotropic tide predicted by the simple 
theory of Bell (1975). Jayne and St. Laurent (2001) found that including this extra term 
in their implementation of the linear SWE (with k ~ 27r/10km) significantly improved fit 
of modeled tidal elevations to those estimated from T/R 


We implemented and tested the scheme of Jayne and St. Laurent (2001), along with two 
variants. The first of these is based on the work of Sjoberg and Stigebrandt [1992], who 
give an alternative (but dimensionally similar) expression for barotropic energy conver- 
sion. In this approach the conversion (comparable to E f of (4)) is calculated by treating 
the bottom topography as a series of discrete steps, and applying theories developed 
for generation of internal waves by flow over sills. Note that in this approach energy 
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conversion is calculated for each step independently, and that this can only be formally 
justified if the steps are far enough apart. Expressions for energy flux away from each 
topographic step are given in Sjoberg and Stigebrandt (1992) and are summarized in 
Gustafson (2001). The resulting expressions for E } again depend quadratically on the 
cross-step tidal velocity, so a linear drag coefficient c/r is readily derived as for the Bell 
formula (4). This second approach has no obvious unknown or tunable parameters, but 
in fact the expression for Ef can be shown to depend strongly on the grid resolution used 
to define the steps (Garrett et al., 2002). For our calculations we used a 1/12° grid to 
define the topographic steps, and stratification profiles from Levitus (1999) to compute 
the IT drag. In numerical experiments with this parameterization it was found necessary 
to allow for an extra tunable scaling factor to obtain the best results. 

The final approach we tried is based on an extension of the theory of Bell to allow for 
two-dimensional topography and a finite-depth ocean (Llewellyn Smith and Young; 2001). 
As with Bell (1975) the theory is linear and inviscid, and small amplitude topography 
is assumed on an otherwise flat bottom. A rigid lid is assumed at the sea-surface. The 
internal tide is forced by the barotropic component of vertical velocity (e.g., Baines, 1982) 
w bt{z) = u • VHz/H where —H < z < 0, and u is the depth independent barotropic 
horizontal velocity. As shown in Llewellyn Smith and Young (2001) the baroclinic bottom 
pressure Pbc(~H) can be given as the convolution of a radially symmetric Green’s func- 
tion Q w (s) and the barotropic vertical velocity at the bottom Pbc(-H) — G cj *w B t(-H). 
Note that in our notation the Green’s function gives the internal wave bottom pressure 
as a function of distance from a “unit magnitude” point source forcing velocity profile (at 
frequency uj) of the form wbt (~) = — z/H. The Green’s function can be written down 
analytically in terms of the flat bottom vertical mode eigenvalues and eigenvectors. A 
slightly simplified version, based on a WKB approximation to the vertical modes can be 
given explicitly as 

g^s) = fe l ~ 1 )Nb Hi (mr(u) 2 - fy^s/HN) , (6) 

n 

where H % is the zero order Hankel function of the second kind, and s denotes radial 
distance (Llewellyn Smith and Young, 2001). The barotropic/baroclinic energy at a fixed 
location is then readily calculated as 

Ef = ( Pbc(-H)w B t(-H )) , ( 7 ) 

where the brackets denote tidal cycle averages. Further details, and exact expressions 
(without the WKB approximation) are given in Llewellyn Smith and Young (2001). 

This theory allows an explicit formulation for the dependence of IT drag on the barotropic 
velocity which is rigorously justifiable for the case of small amplitude bottom topography 
(and laterally homogeneous stratification) 

Ef = ( u • (WH) t Q w (s) * Vifu) = (vl-11 u *u) . (8) 
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Prom (8) we see that for a fixed frequency u> the IT drag component of the dissipative 
stress T in (1) can be represented as ? IT = * U. This is linear in the transports 

V = Hu, but is non-local in space. Note also that this dissipation operator is frequency 
dependent, and hence also non-local in time. As it would be extremely expensive to 
fully implement this in a numerical scheme for solving the time dependent SWE, we 
implemented this scheme only approximately. First, as we focus primarily on the M 2 
tide, we take u> = 1.4052 x 10 4 s 1 as a constant. Second, we do the convolution in space 
once, using the frequency domain barotropic tidal velocity fields u from a numerical model 
to calculate E f (6,(j)) as a function of position using (7). Then, we replace convolution 
with H w by multiplication with the 2 x 2 spatially varying drag tensor 

R(0, 4>) = E f (9, (/>) |u • V# I -2 (VtffVtf. (9) 

Note that R is singular, since conversion to baroclinic motions occurs only for motions 
perpendicular to local bathymetric gradients. 

The approximate linear 2x2 drag tensor R will result in the correct energy dissipation at 
each point provided the tidal velocities in the numerical solution are exactly equal to the a 
priori assumed u in (9). For our calculation we computed u by solving the SWE in a series 
of small overlapping rectangular areas, each 20° on a side, with grid resolution of 1/30°, 
and open boundary conditions from the global inverse solution TPXO.5 (Egbert and 
Erofeeva, 2002). These local solutions are in fact quite similar to TPXO.5. The resulting 
approximate drag tensors R are thus most appropriate for open ocean tidal velocities of 
the modern ocean. Some thoughts on possible refinements of this procedure will be given 
below. The amplitude of the drag tensor (i.e, the sum of the diagonal elements of H~ l R) 
is plotted on a logarithmic scale in Figure 1. The computed IT drag is largest in the open 
ocean over rough topography, in most of the same areas where Egbert and Ray (2000, 
2001) found significant dissipation in the barotropic tide. The IT drag coefficient is also 
often large along the continental shelf break, but note that in these areas tidal volume 
transports, and hence Ef, are typically not so large. 

Note that the spatial pattern of IT drag coefficients are very similar for all three ap- 
proaches tried. This is not surprising since in all cases the estimated energy fluxes scale 
quadratically with topographic slope and tidal velocity, and linearly with bottom strati- 
fication, as in (4). 


3.3 Implementation of the SAL Correction 

In (1) the effects of ocean self-attraction and loading (SAL) are included as an extra 
equilibrium-like tide Csal- This can be related to the tidal elevation £ through convolution 
(on the sphere) with the SAL Green’s function (e.g., Hendershott, 1972; Ray, 1998) 

Csal = Q SAL * C- (10) 
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Substituting (10) into (1) results in an integro-differential equation. A naive solution 
approach, applying the full convolution operator of (10) at each time step, would obviously 
not be computationally practical. For the modern day ocean C is well known (at least 
in the open ocean) so Csal can be computed once and added to (1) as an extra forcing 
term This approach, with elevations from Schwiderski (1978), was used by Le Provost 
et al. (1994). However, solutions computed in this way are not really independent of all 
data, and the computed elevation fields will not in general be consistent with the assumed 
Csal- This inconsistency can lead to significant imbalances in the energy equation (e.g., 
Le Provost and Lyard, 1997). This simple approach is of course even harder to justify for 
calculations with different ocean geometries, where tidal elevations may be expected to 
deviate from the modern estimates of C used to compute Csal- 

Another simple approach that has frequently been used is to approximate convolution 
with Qsal with multiplication by a scalar factor 0. In this case C — Csal is replaced 
by (1 — P)C and (1) is reduced to a partial differential equation. Analysis by Accad and 
Pekeris (1978) suggested 0 « 0.085; Schwiderski (1978) used 0 = 0.1. However, as pointed 
out by Ray (1998) no fixed scalar 0 is appropriate for all locations in the ocean. Our 
own experiments suggested that this scalar approach could result in significant errors in 
global tidal solutions. To allow rigorous treatment of SAL in global modeling Hendershott 
(1972) suggested an iterative approach, with elevations C n from iteration n being used 
in (10) to compute C sali an d then the result used in (1) to compute a new solution 
C n+1 . However, numerical experiments with this approach suggested that convergence of 
this iterative scheme could not be guaranteed (Hendershott, 1972; 1977). We tried this 
iterative scheme for modern ocean basins, starting from an initial guess at C^al estimated 
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from an accurate T/P based model (TPXO.5). In Figure 2a we plot (dashed line) the 
RMS change in elevations |C” — C n-1 | between successive iterates. Differences increase 
with subsequent iterations. Clearly the scheme does not converge in this case. 

A slightly modified iterative approach (which was in fact first suggested by Accad and 
Pekeris (1978)) converges rapidly. The idea is to write 

Csal = PC + ( Qsal * C — PC) = PC + Q'sal * Ci (11) 

where p is chosen to be a reasonable scalar SAL approximation (we used (3 = 0.1). At 
iteration n + lwe then replace C - Csal in (1) with (1 - p)( n+1 - Q' SAL * C n - It is readily 
verified that if this converges as n ->■ oo, C n converges to the solution of the full integro- 
differential equation. The change between successive iterates for this modified scheme, 
which converges in 4-5 iterations, are plotted as the solid line in Figure 2a. Note that 
the first iterate for this modified scheme (C 1 ) was forced with the full Csal computed from 
TPXO.5, and P - 0, so the first iteration is the same as for the simpler scheme. 

In Figure 2b we plot the global integral of work done by the SAL term in the tidal 
equations 


Wsal = {pg J f Csal&C /dt)dS^ . (12) 

Since we assume an elastic Earth (with Love numbers strictly real) Qsal is also real and 
the cycle average of the global integral must be zero (e.g., Egbert and Ray, 2001). W S al 
(computed with C obtained from iteration n and Csal from the previous) is plotted as 
a function of iteration number in Figure 2b. Note that for the initial iteration W S al 
is nearly 0.5 TW, a significant fraction of the 2.5 TW of M 2 tidal dissipation. For the 
modified (but not the simple) iterative scheme the Wsal converges rapidly to zero, further 
demonstrating the consistency of Csal and C obtained in this way. 

The iterative scheme for the SAL convolution suggests a possible refinement of our internal 
tide parameterization. Given a solution at iteration n, the currents u n could be convolved 
with the Green’s functions Q u to compute IT stress for each modeled tidal constituent. 
These in turn could be added as extra forcing terms for iteration n + 1. A modification, 
comparable to that used for the SAL convolution, would be include the local linear drag 
tensor R in the model equations, and only add the deviation from this as a forcing. 
This proposed scheme would (if it converges!) allow for a frequency dependent, spatially 
non-local parameterization for IT drag. Since iteration for the SAL correction is already 
required, this may well require no additional runs of the SWE solver. However, we have 
not tested this scheme. 
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Figure 2: Convergence of iterative schemes for ocean self-attraction and tidal loading, (a) 
RMS difference in elevation between successive iterates, averaged over the globe. Dashed 
line: iterative scheme proposed by Hendershott (1972). Solid line: modified scheme of 
Accad and Pekeris (1978) described in text, (b) Convergence of global integral of SAL 
work term (12) for the simple (dashed lines) and modified (solid lines) iterative schemes. 
For all runs the IT drag parameterization based on (9) was used. 

4 Results 


4.1 Effect of Resolution 

To assess how well the hydrodynamic solutions reproduce modern tidal elevations we 
compare model outputs with the global inverse solution TPXO.5. This reference model is 
most accurate in the open ocean, and equatorward of 66° where T/P data is available. The 
RMS difference between M 2 elevations computed using a wide range of grid resolutions 
and the T/P reference solution are plotted in Figure 3. Increasing resolution of the finite 
difference grid significantly improves agreement between model outputs and elevations 
inferred from T/P. This it true for model runs with and without an internal tide drag 
parameterization, and in deep and shallow water (Figure 3a). Including internal tide 
drag significantly reduces misfits for all resolutions. The best results, with misfits in deep 
water (H > 1000m) slightly below 5 cm RMS, are obtained with internal tide drag r un 
at 1/12 resolution. Only slightly worse results are obtained with a grid resolution of 
1/8 . Including shallow seas in the comparison increases misfits somewhat (thiner lines 
in Figure 3a). The work done by the tidal potential, which must equal the global tidal 
energy dissipation, also converges to approximately the proper value (of approximately 2.5 
TW) as resolution is increased. The agreement is better (especially at lower resolution) 
when the IT drag parameterization is used. 

For the results presented in Figure 3 internal tide drag was parameterized in terms of 
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Figure 3: (a) RMS difference between M 2 elevations from the T/P global inverse solutions 
TPXO.5, and hydrodynamic solutions computed at a range of grid resolutions. Solid and 
dashed lines are for computations with and without internal tide drag respectively. Heavy 
lines marked with circles and pluses give misfits for deep ( H > 1000m) water; the light 
lines give misfits for the full globe (equatorward of 66°). (b) Global integral of work 
done by the M 2 tidal potential, for solutions with (solid lines) and without (dashed lines) 
internal tide drag, for solutions at a range of resolutions. 


the linear drag tensor R defined in (9), without adjusting any parameters. Very similar 
results were obtained with the other two IT drag parameterizations discussed above, but 
in these cases some tuning of the overall scale of the IT drag was required. For the scheme 
based on Bell’s (1972) formula (4) we found results were best for a value of k « 10 -3 m -1 , 
slightly larger than the optimal value k,^2ttx 10 _4 m _1 found by Jayne and St. Laurent 
(2001). For the third scheme IT drag coefficients derived from Sjoberg and Stigebrandt 

(1992) (computed using bathymetry on a 1/12° grid), worked best when multiplied by a 
factor of 0.75. 


4.2 Sensitivity to Errors in Bathymetry 


Figure 3 suggests that further increases in grid resolution will lead to little improvement 
in fit to the observed modern tidal elevations. A number of factors probably limit the 
accuracy that can be expected from numerical solution of the SWE at any resolution. 
For example, the IT drag parameterizations we have tested are all based on a linear 
treatment of small amplitude topography. Even for this simplified linear case a proper 
treatment of IT drag would involve convolution with a non-local operator, so we have 
only incorporated IT drag approximately. The parameterization of bottom drag in terms 
of (3) with a constant c D , or of subgrid scale motions in terms of an eddy diffusivity, must 
also be only approximately correct. And of course, even for the modern ocean, there are 
potentially significant errors in the available bathymetric databases. 
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Case 

Global RMS (cm) 

Deep Ocean RMS (cm) 

5% errors 

3.23 

2.69 

10% errors 

8.17 

6.65 

Variable errors 

8.23 

6.67 


Table 1: Root mean square changes in Af 2 tidal elevations due to random variations of 
bathymetry. 


To test the sensitivity of model solutions to bathymetric errors we did a series of r uns 
with small random perturbations to our standard bathymetry. To save on computer time 
these runs were all done on a 1/4 grid. Several scenarios for the bathymetry errors were 
considered. In each case 10 perturbed bathymetric grids were generated, solutions were 
calculated for each, and the global RMS elevation differences (relative to the standard 
bathymetry solution) were computed. Results are summarized in Table 1. For the first 
case the random variations were 5% of the local depth, with a decorrelation length scale 
of 2.5 . For the second case error amplitudes were increased to 10%, and for the third 
variable error magnitudes were assumed, with a dependence on depth modeled on the 
statistics of differences between the Gtopo30 (Smith and Sandwell, 1997) and ET0P05 
(National Geophysical Data Center, 1992) databases, which were found to be: H < 100m: 
25%; 100 m < H < 200 m: 15%; 200 m < H < 1000 m: 10%; 1000 m< H < 3000 m: 6%; 
H > 3000 m: 3%. °For cases 2 and 3 the decorrelation length scale of the bathymetric 
errors remained 2.5 . 

As another test of sensitivity of model elevations to bathymetric inputs, we ran the model 
with bathymetry derived directly from the ETOP05 database alone. RMS differences of 
M 2 elevations between this solution and that from our standard bathymetry case (again 
computed on a 1/4° grid) were 11.8 cm (all depths), and 6.9 cm ( H > 1000m). This last 
result is at least roughly consistent with the sensitivity inferred from our experiments with 
random depth errors. It is difficult to assess the accuracy of the Gtopo30 bathymetry. 
However, even if errors are only half the size of the difference between Gtopo30 and 
ETOP05 (perhaps an optimistic assessment), Table 1 suggests that open ocean errors 
in modeled M 2 tidal elevations of 3-4 cm RMS should be expected. Given the other 
shortcomings in the numerical SWE model, the 5 cm RMS error achieved by our 1/12° 
model is probably already approaching the lower limit of M 2 tidal modeling errors that 
can be achieved at present. The high sensitivity of model outputs to bathymetry also 
implies that it will be very difficult to model ocean tides accurately in the distant past. 


4.3 M 2 Tides in the Last Glacial Maximum (LGM) 

With the hydrodynamic model tuned to model the present ocean tides, we consider the 
effect of dropping sea— level to that inferred for the LGM. For these experiments Bathy- 
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metric grids were modified using the 1 topography and ice model ICE-4G (Peltier, 1993; 
1994). For each time considered (5, 10, 15, 20 ka), we computed the difference in ocean 
depth with the present day (0 ka) ICE-4G topography. The difference fields were then 
interpolated onto our 1/8 nearly global grid, and added to the present day bathymetry. 
In this manner we made sea-level adjustments consistent with the resolution of ICE-4G, 
but retained the higher resolution bathymetric details of the modern topography. We 
did these numerical experiments at 1/8 (instead of 1/12 ) to economize on computer 
time. In Figure 4a we plot RMS difference between the present day tidal elevations (as 
determined by T/P), and the elevations computed from a series of times over the past 
20 kyr. This figure reveals that M 2 tidal elevations in the LGM were indeed significantly 
different. ° 

In Figure 4b we plot M 2 tidal dissipation (inferred from the astronomical work) in the 
numerical solutions for the past 20 kyr. A significant increase in dissipation, by approxi- 
mately 40% or more, accompanies the drop in sea-level which exposed many of the shallow 
shelf areas where present day dissipation is greatest. 

For the results of Figure 4 we have kept the IT drag tensor exactly as estimated for the 
modern ocean. Since this depends strongly on stratification this extra frictional term 
may not be appropriate for the LGM. (The IT drag based on (9) also depends, but much 
more weakly, on the tidal currents. Note that with the other two formulations tried the 
estimated coefficients depend only on the bathymetric gradient and the stratification.) 
As a sensitivity test we did computations for 20 ka with IT drag coefficients reduced and 
increased (uniformly across the globe) from our estimates based on the modern stratifica- 
tion. The first case, which corresponds to a significant decrease in stratification, results in 
further increases in tidal dissipation to over 4 TW. For the second case (with increased IT 
friction) total dissipation is reduced, and elevations are somewhat more consistent with 
those of the present day. 

In Figure 5 we plot M 2 elevation amplitude and phase for the modern day and 20 ka 
1/8 M 2 tidal solutions, computed with the standard IT drag parameterization, and the 
20 ka solution computed with IT drag coefficients reduced by a factor of 4. Note that 
the color scale for the amplitudes is logarithmic. Amplitudes for the 20 ka solutions 
are noticeably larger throughout the ocean. Amplitude increases are especially great in 
the North Atlantic, where tidal amplitudes in the 20 ka model solution exceed 3-4 m 
over much of the Labrador Sea, and off the Atlantic coasts of Spain and North Africa. 
Reducing the IT drag coefficient results in further increases in amplitude. For this case 
more significant amplitude increases also occur in the Pacific basin, especially off the east 
coast of Australia and around New Zealand. 

In Figure 6 the distribution of tidal energy dissipation is plotted for the three numerical 
solutions of Figure 5, and for the solution computed with present day topography but no 
internal tide drag. For these calculations we followed the approach described in Egbert 
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Figure 4: (a) RMS misfit to modern tidal elevations for solutions computed with 

bathymetry appropriate to the past 20 kyr. (b) Work done by the tidal potential over this 
time. The asterisks and open circles at 20 ka are for the cases with IT drag multiplied by 
0.25 and 2, respectively. 

0 ka 0 ka 20 ka 20 ka 20 ka 

TPXO.5 GOT99hf IT no IT IT ITx.25 ITx2 

Shallow 1.625 1.717 2.104 2.644 2.443 3.088 1.877 

Deep .813 .729 .576 .061 1.513 1.038 1.817 

Total 2.438 2.446 2.680 2.705 3.956 4126 3.694 

Table 2: Dissipation (in TW) in deep and shallow seas, following the division given in 
Egbert and Ray (2001). TPXO.5 and GOT99hf are constrained by T/P data, the other 
4 are purely numerical solutions discussed in text. 


and Ray (2001), with dissipation computed as a local balance between energy flux diver- 
gence and work done by the tide generating force and SAL. A rough breakdown between 
dissipation in shallow seas and the deep ocean, is given in Table 1 for these four solutions, 
and for the 20 ka solution with IT drag increased by a factor of 2. For comparison, results 
from two of the T/P constrained solutions from Egbert and Ray are also given here. A 
more detailed breakdown for the shallow seas and deep-ocean areas discussed in Egbert 
and Ray (2001) is given in Figure 7 for all 5 of the numerical solutions and TPXO.5. 

Figure 7 reveals fairly good agreement between the spatial distribution of dissipation in 
the 0 ka solution (with IT drag) and in the T/P constrained solution TPXO.5. The most 
significant discrepancies in shallow seas are in the Arctic, along the West coast of North 
America and around New Zealand. In all of these areas the numerical model dissipates 
more energy than the T/P based estimates suggest. Note that areas around Indonesia 
are poorly constrained by the T/P data, as discussed in Egbert and Ray (2001). In deep 
water, areas in the South Pacific (Micronesia/Melanesia and Polynesia) come out high in 
the numerical solution, while the Mid-Atlantic and W. Indian Ridges come out low. This 
suggests that our IT parameterization overestimates drag produced over the larger scale 
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Figure 5: Amplitude and phase of M 2 tidal solutions for the present (top) and 20 ka 
(middle), computed with the standard IT drag, and for 20 ka computed with IT drag 
coefficients multiplied by 0.25 (bottom). Note that a logarithmic scale is used for ampli- 
tudes. 
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Figure 6: Spatial distribution of dissipation computed for the numerical solutions as the 
balance between energy flux divergence and work terms. Solution were computed using 
present day bathymetry with (a) and without (b) IT drag parameterization; and for 20 
ka bathymetry with IT drag from the present day (c) and reduced by a factor of 4 (d). 


[§] 


topography associated with volcanic arcs and island chains, and underestimates drag over 
oceanic spreading centers, which are dominated by smaller spatial scales, and are spread 
over a larger area. Turning off the IT drag of course nearly eliminates dissipation in all 
of the deep ocean areas. It also causes significant increases in dissipation in a number of 
shallow seas, especially in area (1), which includes the Labrador sea and all points north 
and west into Hudson Bay, on the Patagonian and European Shelves, and in the Bering 
Sea. 

As Figures 6 and 7 and Table 2 show, dramatic changes to the distribution of tidal 
energy dissipation are likely to result from the drop in sea-level of the LGM. Some of 
the major shallow sea sinks from the present day ocean are significantly reduced in area, 
and as a result would dissipate little energy. These include the Yellow Sea, the shelf 
off the northeast coast of Brazil, and the Andaman Sea. Dissipation is also significantly 
reduced over the European and Patagonian shelves and in the Gulf of Maine off the east 
coast of North America. Dissipation in some other areas increases dramatically in the 
model solutions for 20 ka. These include the Hudson Bay/Labrador Sea area, the Arctic 
(including here the Norwegian sea) and Antarctic, and the area around the Caribbean. 
The dissipation in this last area is concentrated over the volcanic arc on the eastern edge 
of this area, and almost all of the dissipation is in fact associated with IT drag, rather 
than bottom boundary layer processes in shallow seas. This is just one of the areas around 
the North Atlantic where tidal dissipation increases dramatically in the LGM. There are 
also significant increases for the Mid-Atlantic Ridge (especially in the North Atlantic; see 
Figure 6), and off the northwest coast of Africa. 

When IT drag is also reduced by a factor of 4 for the 20 ka bathymetry, further significant 
increases are observed in the Hudson Bay region, the Arctic and Antarctic, and on the 
Eastern edge of the Caribbean. Each of these areas now accounts for a rather astounding 
(by the standards of the present day ocean) 400-600 GW. Dissipation is reduced in some, 
but not all, deep ocean areas. Thus significant reductions in the drag coefficients are 
counterbalanced to some extent by increases in tidal current velocities. Thus while the 
IT drag coefficient is reduced by a factor of 4, dissipation in the deep ocean areas (where 
IT drag is by far dominant) is only reduced by a factor of 1.5 overall (Table 2). Figure 6d 
provides a striking image of the redistribution in tidal dissipation that would result from 
this (rather hypothetical) change in IT drag coupled with the drop in sea-level. In this 
case there are three centers of tidal dissipation, in the North Atlantic, in the seas around 
New Zealand, and around Antarctica. 

Increasing IT drag by a factor of 2 for the 20 ka bathymetry results in less obvious changes 
in the distribution of dissipation. However dissipation in the deep water is increased 
somewhat (Table 2), especially along the Mid- Atlantic and West Indian Ridges. 
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Figure 7: (a) Energy dissipation integrated over shallow seas and selected deep ocean areas 
defined in (b). 
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5 Conclusions 


Efforts to develop a purely hydrodynamic model of the present day tides have been rea- 
sonably successful. A convergent iterative scheme for rigorous treatment of SAL has been 
developed, and three different approaches to parameterization of IT drag have been tried. 
Both the rigorous treatment of SAL and IT drag parameterizations are required to ob- 
tain agreement with the highly accurate T/P constrained tidal solutions for the present 
day ocean. Success with the numerical modeling also required high resolution numerical 
grids and accurate bathymetry. Limitations in our ability to reproduce the present day 
tides are probably due to errors in bathymetry, and probably also our approximate linear 
theory for computing IT drag. 

After achieving reasonable accuracy for the present day tides, the hydrodynamic model 
was run with bathymetry estimated for a series of times over the past 20 ka, extending into 
the LGM. The drop in sea-level associated with the LGM is found to result in significant 
changes in tidal fields, and increases of 40-50% in total dissipation. However, details are 
sensitive to stratification, and this is much more poorly known than bathymetry in the 
recent past. The full implications of the predicted changes in LGM tides remain to be 
explored. 
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